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The  solution  to  a  differential  equation  of  the  form 

x  =  Ax(t),  a;(0)  =  xq 

is  the  function  x(t)  =  eAtXo  [5].  The  expression  eAt  is  the  matrix  exponential  function.  Examples  of  such  equations 
arise  in  control  theory  and  tracking  applications. 

A  key  application  is  the  tracking  of  a  ballistic  target  using  noisy  measurements.  In  this  case,  the  matrix  A  above 
is  actually  a  non-linear  function  of  both  x  and  t.  The  extended  Kalman  filter  (EKF)  has  been  used  in  these  tracking 
applications  [1,  2],  The  typical  formulation  of  the  EKF  uses  a  first  or  second-order  approximation  to  the  solution  of 
the  differential  equation  to  save  operations  [3],  While  such  implementation  is  efficient,  it  has  been  shown  that  in  some 
conditions  the  EKF  may  show  significant  bias  in  altitude  and  ballistic  coefficient  [6],  Under  such  conditions  it  may  be 
preferable  to  use  the  matrix  exponential  function  directly. 

In  this  paper  we  describe  and  benchmark  an  implementation  of  the  matrix  exponential  function.  The  implemen¬ 
tation  is  based  on  the  standard  technique  of  “scaling  and  squaring”  from  the  literature  [4,  5],  The  major  kernels  in 
this  technique  are  matrix  multiplication  and  Gaussian  elimination.  In  the  matrix  multiply  kernel,  the  implementation 
makes  use  of  SIMD  vector  extensions  present  on  the  PowerPC  G4  (Altivec)  and  the  Intel  Xeon  (SSE-2).  Although  the 
use  of  the  matrix  exponential  expands  the  operation  count  of  the  extended  Kalman  filter  substantially,  benchmarks  of 
the  implementation  show  that  the  workload  is  well  within  the  capabilities  of  modern  processors. 
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Definition 

[Moler  and  Van  Loan, 2003] 


The  solution  to  the  differential  equation 

x  =  Ax(t) 

x(0)  =  x0 


is  given  by 

x(t)=  eAtx0 

At 

Where  e  is  the  matrix  exponential  function, 

At  r  .  A2t2 

C  —  /+  At  H - h  ... 

2! 


Notice  that  if  A  =  [a,;],e''  *  [ev]in  general. 
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Application:  Ballistic  Target  Tracking 


*  Tracking  of  a  ballistic  target  using  noisy  measurements 

*  Tracking  accomplished  using  the  extended  Kalman  filter 

-  “extended”  means  that  system  dynamics  are  non-linear 
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The  Extended  Kalman  Filter 


E 


ill 


Estimate  next  state  based  on  previous  state  and  new  measurement 


Measurement 


Updated 
~+  state  and 
xk  ’  *k  covariance 
estimates 


K  =  +  <2(0, 

where  O  =eJ'  for  a  matrix  J, 

and  Q(t)  is  the  process  noise  covariance. 
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Calculation  Overview 
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Preferred  method,  Pade  approximation,  is 
only  valid  when  IIAII  is  small 

Use  the  fact  that  eA  =  ( eAlm)m 


1 .  Choose  an  integer  j  and  scale  A  by  m= 

2.  Use  a  Pade  approximation  to  calculate  E  =  eAiV 

3.  Perform  j  matrix  multiplies  to  calculate  E~ 


This  technique  is  referred  to  as  “scaling  and  squaring”  [4,5]. 
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X  =  A; 
c  =  1; 

E  =  I; 

D  =  I; 

for(k  =  1;  k  <=  q; 

{ 

c  =  c  *  (q-k+1) 
X  =  A*X; 

E  =  E  +  cX; 
if  (k  is  even) 

D  =  D  +  cX; 
else 


D  =  D  -  cX; 

} 

E  =  D\E; 


k++)  //  q=number  of  iterations 

/  (k* (2*q— k+1) ) ; 

//  Matrix  multiply 
//  Matrix  scale  and  add 
//  Matrix  add  or  subtract 


//  Solve  using  LU  factorization 
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Implementation  Overview 


Step 

Operations 

Percentage 
of  op  count 

Scale  the 
matrix  A 

Elementwise 

multiply 

<2% 

Pade  iteration 

Matrix  multiply, 
scale,  add 

50-75% 

LU  and 
backsolve 

3-6% 

Repeated 

squaring 

Matrix  multiply 

13-50% 

Op  counts  assume  6  Pade  iterations 


Implementation  Features 

*  Single-precision  real  or 
complex  float 

*  C++ 

*  Uses  an  object  for 
storage 

*  Calls  VSIPL  routines 

*  Uses  Altivec-optimized 
matrix  multiply 

*  Choose  accuracy  to 
match  limits  of  single¬ 
precision  calculations 


void  create (Matrix<T>  &A, 

//  Allocates  memory  &  initializes 

Matrix<T>  &E) ; 

//  LU  factorization 

void  run (  Matrix<T>  &A, 

//  Performs  computation 

Matrix<T>  &E)  ; 
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Performance 


Achieved  Performance  on  PowerPC  G4 


# 

> 

U 

c 

© 


Complex 

Real 


*  Platform:  Mercury  500  MHz  PowerPC  G4 

*  Achieves  respectable  performance  for  large  matrices 

*  For  tracking,  sizes  of  interest  are  small  —  6x6  matrices 

-  A  tuned  implementation  could  be  produced  for  this  size 

MIT  Lincoln  Laboratory 
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Performance  Breakdown 


Matrix  size 


■  Squaring 

□  LU  and  backsolve 

■  Pade  iteration 

□  Scale 


Matrix  size 

Scale  — A — Pade  iteration  LU  and  backsolve  — * — Squaring 


Performance  breakdown 
on  PowerPC  G4 

Steps  based  on  matrix 
multiply  are  more 
efficient  than  other  steps 

For  large  matrices,  matrix 
multiply  steps  still 
consume  most  of  the 
execution  time 

LU/backsolve  is  a 
substantial  percentage  of 
time  despite  being  a  low 
percentage  of  the  op 
count 
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The  Matrix  Exponential  in  Tracking 
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Trackable  targets 


Matrix  exponential  is  a 
substantial  part  of  the  EKF’s 
operation  count 

How  many  targets  could  a 
single  processor  track? 

-  Assume  500  MHz  PPC  G4 

-  Use  execution  time  of  6x6  real 
matrix  exponential 

-  Assume  remainder  of  EKF 
has  efficiency  comparable  to 
LU  factorization  (~0.04%) 

-  Vary  track  rate  from  2-10  Hz 

A  single  processor  can 
potentially  track  many  targets 


MatrixExponential-1 0 
JML  9/20/04 


MIT  Lincoln  Laboratory 


Conclusions 
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*  Matrix  exponential  function  is  important  for 
tracking  applications 

*  A  large  percentage  of  the  operations  are  matrix 
multiply  functions 

*  An  efficient  implementation  of  this  function  allows 
it  to  be  used  in  an  extended  Kalman  filter 

*  Many  targets  can  be  tracked  using  even  a  single 
processor 

-  Using  multiple  processors  obviously  allows  more 
targets  to  be  tracked 
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Application:  Ballistic  Target  Tracking 


*  Tracking  of  a  ballistic  target  using  noisy  measurements 

*  Tracking  accomplished  using  the  extended  Kalman  filter 

-  “extended”  means  that  system  dynamics  are  non-linear 
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The  Matrix  Exponential  in  Tracking 
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Trackable  targets 


Matrix  exponential  is  a 
substantial  part  of  the  EKF’s 
operation  count 

How  many  targets  could  a 
single  processor  track? 

-  Assume  500  MHz  PPC  G4 

-  Use  execution  time  of  6x6  real 
matrix  exponential 

-  Assume  remainder  of  EKF 
has  efficiency  comparable  to 
LU  factorization  (~0.04%) 

-  Vary  track  rate  from  2-10  Hz 

A  single  processor  can 
potentially  track  many  targets 
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